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1.  Introduction 

The  recent  success  of  characteristic  based  upwind  differencing  on  structured  meshes  has  spurred 
significant  interest  into  adopting  such  methods  into  unstructured  adaptive  mesh  aigo- 
rithmsf  Both  adaptation  and  upwind  differencing  strive  to  increase  the  accuracy  and 

efficiency  of  the  numerical  simulation,  and  the  combination  of  the  two  approaches  suggests  a 
powerful  and  flexible  platform  for  the  simulation  of  flows  with  features  spanning  a  wide 
spectrum  of  length  scales. 

Adaptive  mesh  methods  and  upwind  techniques  were  both  developed  and  tested  on  inviscid  flow 
problems,  and  it  may  be  convincingly  argued  that  some  aspects  of  their  application  to  multi¬ 
dimensional  viscous  flow  remain  incompletely  understood.  Most  successful  upwind  techniques, 
for  example,  rely  upon  the  solution  of  a  Riemann  problem  on  a  direction-by-direction  basis^l. 
The  spatial  operator  resulting  from  this  approach  cannot  be  considered  truly  isotropic  and  its 
interaction  with  the  viscous  discretization  (often  central  difference)  has  not  been  extensively 
researched.  Moreover,  in  viscous  dominated  regions,  the  Riemann  problem  solved  by  the 
convective  discretization  becomes  ill-posed  due  to  the  presence  of  large  viscous  terms  in  the 
governing  equations.  Nevertheless,  the  overwhelming  success  of  upwind  methods  at  providing 
high  quality  solutions  to  complex  flows  continues  to  motivate  considerable  research  into  these 
and  other  fundamental  topics^5!. 
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Unstructured,  adaptive,  mesh  methods  have  enjoyed  a  similar  degree  of  success  in  the  last 
decade.  From  their  initial  application  to  inviscid  2D  fiowsl6l-!7l.l8],  through  complete  three 
dimensional  simulations  over  complex  geometriest9!-!  101.11  ll(  such  techniques  have  promised 
outstanding  efficiency  by  clustering  mesh  nodes  where  they  are  most  needed  in  the  evolving 
solution.  Nevertheless,  several  outstanding  issues  still  exist  Credible  and  general  techniques 
for  deciding  where  to  adapt  the  solution  have  yet  to  be  definitively  determined,  and  fundamental 
questions  exist  concerning  mesh  convergence  and  consistency!12!.  Current  research  efforts  are 
also  focusing  ever  more  sharply  on  issues  of  mesh  quality  and  the  local  disruption  of  mesh 
metrics  caused  by  introducing  new  nodes  into  the  computational  domain!  Dl.f  14]  The  issue  of 
mesh  quality  is  of  particular  importance  in  the  discretization  of  the  viscous  terms  in  the  full 
Navier-Stokes  equations.  Hie  second  differences  in  such  terms  are  computed  as  the  first  differ¬ 
ence  of  first  differences,  thus  compounding  mesh  stretching  or  skewing  errors  and  resulting  in 
much  more  stringent  requirements  on  mesh  or  element  quality.  This  problem  is  further  com¬ 
pounded  by  the  highly  nonlinear  nature  of  viscous  flow  which  increases  the  importance  of 
resolving  subtle  flow  features  which  may  be  very  readily  disturbed  by  mesh  irregularity.  As  a 
result,  conclusive  demonstration  of  adaptive  unstructured  mesh  techniques  in  resolving  fine 
detail  within  an  adapted  flow  field  remains  elusive. 

Recently,  an  upwind  based  solver  for  the  3D  Navier-Stokes  equations  was  introduced!*5!.  The 
method  uses  adaptively  refined  hexahedral  based  meshes  and  permits  adaptation  through 
directional  cell  division.  Such  cells  provide  a  natural  environment  for  mapping  the  multi¬ 
dimensional  upwind  stencil  which  relies  on  ID  operators  in  each  spatial  direction.  Additionally, 
the  tessellation  easily  forms  the  high  quality  surface  mesh  required  for  accurate  evaluation  of  the 
stress  tensor  and  heat  flux  vector. 

This  report  focuses  on  the  viscous  aspects  of  the  technique  and  includes  a  detailed  presentation 
of  the  discretization  and  implementation  of  these  terms.  This  includes  a  discussion  of  the  order 
of  accuracy  and  its  degeneration  on  nonuniform  meshes.  Tt  also  includes  an  examination  of 
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several  basic  numerical  properties  designed  into  the  scheme.  A  section  addressing  mesh 
smoothness  and  surface  point  introduction  is  also  included.  Various  sub-  and  transonic  example 
problems  facilitate  a  discussion  of  the  topics  mentioned  in  the  preceding  paragraphs.  These 
examples  also  provide  a  basis  for  examining  the  interaction  of  the  upwind  scheme  with  the 
centrally  differenced  viscous  terms  -  especially  within  the  boundary  layer  itself  and  near  regions 
where  the  entropy  conditional  may  be  violated. 
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2.  Description  of  Method 


2.1  Governing  Equations 

The  Navier-Stokes  equations  describing  the  unsteady  flow  of  a  viscous  perfect  gas  may  be 


written  in  integral  form: 


F&F[~Fv 


U) 

{2} 


Here  V  is  the  state  vector  of  conserved  variables,  and  F  is  the  complete  tensor  of  flux  density 


which  contains  both  inviscid  and  viscous  components.  V  refers  to  an  arbitrary  control  volume 
and  dV  is  its  closed  boundary  with  the  outward  facing  unit  vector  n  =  [nx^iynz]^.  In  Eqs.{  1 } 
and  {2}: 


p 

pul 

+ 

PV  J 

+ 

pw  k 

pu 

(pu2  +p)i 

+ 

P«V  j 

+ 
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Oi  +  Oj  +  Ok 

Xxx\  +  rxyj  +  xxzk 

Tyxi  +  tyyi  +  Xy  Zk 

?zxi  +  ?zy  j  +  \ zr ^ 

nxi  +  Fiyj  +  nz  k 


nx=UTxx  +  vTxy  +  H>T«  - 

Tly^tXyx  +  Vtyy  +  W  Xyl  ~  Ijy 

where 

FFz-uXzx  +  vXzy  +  wxzz~qz 


In  this  work,  the  equation  of  state  for  a  calorically  perfect  gas  relates  the  mechanical  and 
thermodynamic  properties  of  the  fluid 

p=(r-i)(p£-eli)  141 

The  viscous  part  of  the  flux  density  tensor  makes  use  of  the  shear  stress  tensor  and  heat  flux 
vector.  These  are  modeled  by  adopting  the  Stokes  Hypothesis,  while  assuming  an  isotropic, 
Newtonian  fluid  yields  a  symmetric  shear  stress  tensor. 

*xx  =  2lLUx  -  Ux  +  Vy  +  Wz) 

Xyy  ~  2 )lVy  ~  UX  +  Vy  +  Wj ) 

%  =  2 tiwz  -  JK UX  +  Vy  +  wj 
Xxy-  XyX  =  H(Uy  +  Vx) 

*xi  =  *zx  =  *“(«z  +  wx) 

XyZ  =  %=M(VZ  +  Wy) 


qy  ~ 
qz~ 


~k 


-k 


BT 

Bx 

BT 

By 

BT 

Bz 


{5} 


Sutherland’s  law  approximates  the  fluctuation  of  molecular  viscosity  with  temperature. 


/i  = 


!  fr.+  110.4 
Mt  +  no.4. 


(T  in  Kelvin) 


{6} 


Finally,  the  Prandtl  number  relates  the  thermal  conductivity  to  viscosity. 


k  = 


y-  1  Pr 


m 
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2.2  Spatial  Discretization 

The  discretization  of  these  equations  may  be  conveniently  developed  by  focusing  on  the 
symbolic  form  of  the  governing  equations  given  in  expression  { 1 } .  Considering  a  control 
volume  fixed  in  time  and  applying  the  mean  value  theorem,  this  equation  may  be  recast  as: 

(sH  =  -vj  r'ndS  (81 

dV 

Equation  { 8 }  may  be  further  specialized  for  application  to  a  specific  polyhedral  control  volume, 
constructed,  with  planar  faces,  around  node  i  (Fig.  1).  Equation  {8}  may  then  be  approximated 
by: 


where  k  denotes  the  fdh  face  of  V„  ~Sk  =  [s**,  Syk,  szk]T  is  the  surface  vector  of  face  k.  The  first 
term  on  the  RHS  of  Eq.{9}  performs  the  inviscid  surface  integral  while  the  second  term  balances 
the  viscous  fluxes  through  the  faces  of  V,\ 

The  solution  method  separates  the  modeling  of  convective  and  viscous  fluxes.  It  uses  an  upwind 
representation  of  the  inviscid  fluxs,  while  the  viscous  fluxes  utilize  central  differencing.  The 
convective  modeling  has  been  extensively  documented  and  Refs.  [17],  [2],  and  Ref.  [19]  contain 
details  of  its  formulation  and  validation.  Ref.  [20]  contains  structured  mesh  results 
demonstrating  the  exceptional  shock  capturing  properties  of  this  formulation,  while  Ref.  [21] 
demonstrates  that  this  inviscid  discretization  also  accurately  represents  smooth  features  with  a 
low  level  of  background  diffusion  (as  compared  to  both  central  and  other  upwind  methods). 
Ref.  [2]  focuses  solely  on  the  unstructured  implementation  of  this  method  and  shows  that  the 
properties  associated  with  this  scheme  on  a  structured  mesh  are  retained  in  the  unstructured 
procedure.  Additionally,  the  inviscid  discretization  has  a  central  difference  option  using  the 
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Figure  1.  Formation  of  dual  mesh  of  Auxiliary  Cells  in  two  and  three  dimensions. 
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blended  2nd  and  4th  order  artificial  dissipation  from  Ref.  [22].  This  permits  a  direct  comparison 
of  upwind  and  central  results  on  identical  adapted  meshes. 

The  integration  of  both  inviscid  and  viscous  fluxes  proceeds  on  a  dual  mesh  formed  by 
connecting  the  geometric  centers  of  the  cells  of  the  physical  mesh.  Figure  1  illustrates  the 
formation  of  this  dual,  auxiliary  mesh,  in  two  and  three  dimensions.  The  inviscid  and  viscous 
fluxes  are  balanced  through  the  faces  of  the  auxiliary  cells  resulting  in  a  node-based  method 
which  promotes  accurate  and  straightforward  treatment  of  boundaries.  Evaluating  the  summa¬ 
tions  of  Eq.{9}  then  forms  an  update  of  the  state  vector  at  each  node  through  a  multi-stage 
Runge-Kutta  time  integration  scheme. 

2.3  Discretization  Strategy 

As  suggested  by  the  form  of  Eq.{9],  the  complete  update  to  any  node,  A{7„  is  a  summation  of 
inviscid  and  viscous  contributions. 

AUi^ATJh^AUhv  {10} 

The  choice  of  central  differencing  for  the  viscous  fluxes  allows  the  scheme  to  be  written 
extremely  compactly.  Eventually  it  involves  only  node-to-cell  and  cell-to-node  communication 
while  organizing  all  gather-scatter  operations  so  that  they  may  be  easily  grouped  for  rapid 
processing.  These  communication  issues  remain  important  because  the  current  implementation 
uses  a  cell  based  data  structure  and  each  physical  cell  can  address  only  the  nodes  at  its  vertices 
directly  WM23!.  Thus,  an  implementation  based  upon  cell-to-node  operations  is  especially 
convenient 

It  is  worth  noting  that  this  basic  discretization  is  also  amenable  to  edge-based  data  structures.  In 
fact.  Ref.  [24]  describes  the  formulation  of  Hessian  and  Laplacian  operators  using  only  edge- 
based  formulas.  For  the  hexahedral  based  meshes  discussed  here,  however,  this  strategy  holds 
little  obvious  advantage,  and  the  cell-based  formulas  were  used. 


In  two  dimensions,  the  viscous  discretization  degenerates  into  a  form  which  is  conceptually 
similar  to  that  first  presented  by  Ref.  [23].  The  current  work  contains  some  important 
differences  in  the  handling  of  surface  vectors,  time  steps,  and  volumes,  however,  and  the 
extension  to  three  dimensions  is  completely  new. 

Figure  1  contains  a  view  of  the  cells  surrounding  node  i  in  two  and  three  dimensions.  The  stress 
tensor  and  viscous  fluxes  are  computed  on  each  of  the  faces  ( N,S,E,W,F,B )  of  auxiliary  cell  i  to 
form  its  viscous  update  AUvisc. 

AU^sc  =  ~  X 

y  kedV 

Applying  the  mid-point  rule  for  integration  and  suppressing  the  subscript  (y)  for  compactness: 

AUvisc  =  Y^N^-  W'S*  (ii) 

+  FF*Sf-  fB»SB] 

The  superscripts  in  Eq.{  1 1 }  refer  to  the  face  associated  with  each  quantity,  and  thus  is  the 
surface  vector  of  the  North  face  of  the  auxiliary  cell  surrounding  node  i.  The  sign  convention 
used  takes  the  Cartesian  components  of  the  surface  vectors  as  positive  when  oriented  in  positive 
coordinate  directions.  These  surface  vectors  are  also  necessary  for  the  inviscid  integration  and 
are  computed  only  once. 

Each  of  the  viscous  fluxes  in  Eq.  {11}  is  constructed  through  a  linear  combination  of  elements 
of  the  shear  stress  tensor  (see  Eqs.{3}  and  {5}).  Filling  this  tensor  requires  spatial  derivatives 
of  velocity  and  temperature.  These  derivatives  may  be  conveniently  expressed  in  tensor  notation 
by  defining  the  free  index  £  =  1,2,  or  3  and  denoting  the  x,y,z  axes  as  xj,  X2,  and  xj. 
=  d{( )  denotes  partial  differentiation  with  respect  to  the  coordinate  directions.  Taking  a 

second  free  index  tj  =  1,  2,  3,  or  4,  the  vector  $  may  be  defined  as  $  =  [u,v,w,r|T,  any 
element  of  which  may  be  referred  to  as  0jj.  With  these  conventions,  the  complete  set  of  first 
derivatives  needed  to  fill  the  shear  stress  tensor  may  be  denoted  as  simply  d^. 
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du  dv  dw 

dx  dx  dx 

du  dv  dw 

dy  dy  dy 

du  5v 
dz  dz  dz 

These  derivatives  must  be  computed  on  all  the  faces  of  each  auxiliary  cell,  and  they  are  evaluated 
by  taking  surface  integrals  over  the  secondary  cells  which  surround  the  centers  of  the  faces  of 
auxiliary  cells.  Figure  2  illustrates  the  construction  of  one  such  secondary  cell  for  the  East  face 
of  auxiliary  cell  i  in  two  and  three  dimensions. 

Each  face  of  the  auxiliary  cell  is  uniquely  associated  with  one  of  the  edges  incident  upon  node  i. 
Thus,  the  bracket  in  Eq.  {11}  contains  one  term  for  every  edge  which  terminates  at  t,  and  we 
actually  seek  to  evaluate  d ^  on  the  midpoint  of  each  edge  in  the  physical  mesh.  This  point 
becomes  increasingly  important  when  considering  more  general  polyhedral  control  volumes  (i.e. 
at  interfaces  of  adapted  regions,  or  in  meshes  composed  of  nonhexahedral  cells).  This  obser¬ 
vation  make  obvious  the  fact  that  exactly  one  secondary  cell  surrounds  each  edge  of  the  physical 
mesh.  Furthermore,  construction  of  the  surface  vectors  for  the  secondary  cells  is  now  very 
convenient,  since  they  may  be  obtained  through  simple  averaging  of  those  associated  with  the 
auxiliary  cells  at  either  end  of  the  edge.  This  construction  requires  only  one  addition,  and  one 
multiplication  per  face  and  never  more  than  edge-to-node  communication.  Figure  3  illustrates  the 
formation  of  the  East  secondary  ceil  around  edge  y.  The  proof  in  the  Appendix  shows  that  this 
construction  will  always  result  in  closed  polyhedra. 

With  the  secondary  cells  defined.  Gauss’  Theorem  provides  an  expression  for  the  first 
derivatives,  along  all  edges  of  the  physical  mesh. 

«  hv  {13} 


ar 

dx 

dT 

dy 

dT 

dz 


where 


4=  12,or3 
rj=  l,2,3,or4 


{12} 
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Figure  3.  Construction  of  surface  vectors  on  secondary  cell  around  edge  T}. 
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where  is  a  component  of  the  outward  facing  normal  vector  of  unit  length,  and  the  volume  is 
that  of  the  secondary  cell.  For  the  secondary  cell  surrounding  edge  y  in  Figure  2  this  becomes: 

d&n  =  2  vt  +  **  + 

+  +  )  - 0"  (5^  +  S£ )  . ,  4 , 

+  +S&)~  *i*(‘ SSi  +  Sij  )] 

The  superscripts  on  Qrj  refer  to  the  location  at  which  the  scalars  are  evaluated  (from  Fig.  2). 
Similar  expressions  may  be  constructed  on  the  other  edges  incident  on  node  i  which  facilitate  the 
formation  of  the  viscous  fluxes  at  all  faces  of  the  auxiliary  cell  surrounding  node  i  and, 
ultimately,  the  evaluation  of  Eq.  {11}  for  AUvisc  at  i. 

Since  it  is  based  on  standard  second-order  central  differencing,  the  discretization  formally  retains 
this  order  of  accuracy  on  smooth  meshes.  Furthermore,  Ref.  [23]  shows  that,  in  2D,  the 
overlapping  secondary  cells  results  in  a  discretization  which  is  capable  of  recognizing  and 
damping  saw-tooth  oscillations  in  the  solution.  Finally,  the  Appendix  demonstrates  that  the 
current  choice  of  auxiliary  and  secondary  control  volumes  results  in  an  integration  which 
preserves  a  zero  gradient  flow  on  any  arbitrary,  nonoverlapping  mesh  (preserves  free  stream 
flow). 

2.4  Unstructured  Implementation 

An  examination  of  the  basic  viscous  discretization  outlined  in  Eqs.{  1 1 }  and  { 14}  reveals  that  it 
involves  only  the  cells  and  nodes  immediately  adjacent  to  node  i.  Furthermore,  these  are  simply 
central  differences  and  the  whole  procedure  may  be  divided  into  a  sequence  of  steps  -  with  no 
single  step  requiring  more  than  nearest  neighbor  communication.  The  goal  of  this  section  is  to 
isolate  the  contribution  of  a  particular  physical  cell  to  the  viscous  stencils  of  each  of  its  vertices. 
This  will  permit  the  formation  of  the  viscous  updates  at  all  mesh  nodes  by  a  single  sweep 
through  the  physical  mesh  cells. 


Integration  of  the  secondary  cell  surrounding  edge  y  requires  operations  in  physical  cells  C,  D, 
Gr  and  H  to  provide  d^rj  on  the  East  face  of  the  auxiliary  cell  around  node  i.  Figure  4 
highlights  the  portion  of  the  integration  (in  Eq.{  14})  contributed  to  by  cell  C.  The  other 
physical  cells  surrounding  this  edge,  D,  G,  and  H,  contribute  to  this  surface  integral  in  a  similar 
fashion.  Assuming  the  volumes  of  the  physical  cells  are  all  nearly  equal  permits  isolation  of  the 
contribution  of  physical  cell  C  to  the  tensor  of  first  derivatives  dgpjj  along  edge  ij. 

Assume  Vc  «  V©  *  Vfc  *  Vfo  ®  V 
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This  edge  is  the  front-south  (fs)  edge  of  physical  cell  C  and  the  notation  {d^pn  )£  reads  as  “the 
contribution  to  dgf>n  on  the  front-south  edge  of  physical  cell  C.”  Figs.  1  and  2  identify  the 
locations  referred  to  by  the  superscripts  on  while  Fig.  4  helps  to  clarify  the  edge  labels  of  the 
physical  cells.  The  zeros  left  in  this  equation  facilitate  direct  comparison  with  Eq.{  14}  and  the 
approximation  preceding  Eq.{  15}  is  necessary  to  separate  out  the  individual  contribution  from 
the  cells  which  border  edge  y.  This  approximation  becomes  exact  on  uniform  meshes,  and  on 
sufficiently  fine  smooth  meshes,  so  it  does  not  alter  the  formal  second-order  accuracy  of  the 
method.  On  stretched  meshes,  one  may  show  that  the  additional  truncation  error  created  by  this 
approximation  is  0( Ax2)  and  is  smaller  than  the  first  order  error  created  by  the  mesh  stretching 
itself. 


In  the  same  manner  that  cell  C  contributes  to  a  portion  of  the  secondary  cell  integration  of  node  /, 
it  also  contributes  to  a  part  of  the  integration  of  the  viscous  fluxes  on  the  auxiliary  cell  around  i 
(Eq.{  1 1 }).  After  some  algebraic  manipulation,  C’s  role  in  the  flux  balance  of  auxiliary  cell  i 
may  be  isolated.  Since  this  expression  involves  the  complete  viscous  part  of  the  flux  density 
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Figure  4.  Contribution  of  physical  cell  C  to  the  integration  of  the  secondary  cell  surrounding 
edge  Tj  -  the  front-south  (fs)  edge  of  cell  C. 
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tensor  (denoted  by  simply  F  for  compactness),  it  is  most  convenient  to  express  it  using  the 
symbolic  form  of  the  notation  used  previously: 

(tWi  •Sf'+Ffc*Sf-  Fcs  •  S/»]  1 16> 

The  surface  vectors  on  the  right  of  this  expression  refer  to  the  faces  of  the  auxiliary  ceil  around  i 
since  this  is  the  cell  being  integrated,  and  the  superscripts  and  subscripts  on  F  retain  their  use 
from  the  previous  equation.  represents  the  contribution  from  ceil  C  to  the  viscous  part  of  the 
flux  density  tensor  along  cell  C’s  front-south  edge.  The  bracket  on  the  right  of  Eq.{  16} 
contains  three  terms  -  one  for  each  edge  of  cell  C  which  is  incident  upon  node  i.  As  these 
expressions  are  applied  to  nonhexahedral  control  volumes,  this  property  is  retained. 

Eq.{16}  is  really  a  Distribution  Formula  -  in  the  cell-vertex  sense.  That  is,  it  contains  a 
contribution  from  all  operations  within  cell  C  that  apply  to  the  viscous  update  at  node  i  All  these 
operations  are  performed  completely  and  without  duplication  within  ceil  C,  and  they  are  then 
distributed  to  node  i.  When  this  node  accumulates  all  the  distributions  from  its  surrounding 
cells,  tne  complete  viscous  discretization  stencil  forms.  Moreover,  information  requests  are  re¬ 
stricted  to  celi-to-node  or  node-to-cell  inquiries. 

In  a  similar  manner,  all  the  physical  cells  contribute  to  each  of  their  vertices.  The  complete  set  of 
distribution  formulae  for  a  hexahedral  cell  such  as  C  are: 


(Wvisc  =  %  *  rf  -  F?  •  S? -  F?  •  SP)  {17} 
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When  these  expressions  are  evaluated  in  all  the  physical  cells  of  the  domain,  the  complete 
viscous  update  is  formed  at  all  the  nodes. 

2.5  Eigenvalue  Scaling  for  High  Aspect  Ratio  Cells 

Near  wall  boundaries,  high  aspect  ratio  cells  are  commonly  used  to  efficiently  resolve  the 
boundary  layet.  Since  the  artificial  dissipation  is  scaled  isotropically  with  the  spectral  radii  of  the 
Jacobian  matrices,  the  damping  in  the  wall  normal  direction  may  become  excessive.  Ref.  [25] 
proposed  a  3D  extension  to  the  2D  variable  scaling  of  Martinelli  [26J  which  adjusts  the  levels  of 
the  blended  2nd  and  4th  order  smoothing  to  be  used  in  conjunction  with  the  central  difference 
convective  fluxes.  When  the  TVD  convective  discretization  is  chosen,  this  non-isotropic  scaling 
affects  only  the  waves  which  are  subjected  to  the  entropy  cutofff  jn  such  cases,  the  scaling 
takes  a  form  similar  to  that  found  in  Ref.  [27]. 
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3.  Adaptation 

Adaptation  increases  the  accuracy  of  the  discrete  solution  by  locally  reducing  the  local  mesh 
dimension  through  directional  division  of  mesh  cells  (physical  cells).  This  process  has  evolved 
into  a  two-stage  procedure.  The  detection/division  stage  begins  by  scanning  a  preliminary 
discrete  solution  for  regions  of  interest.  It  then  enhances  the  mesh  in  those  regions  through  cell 
division  to  create  new  computational  nodes.  The  second  stage  consists  of  mesh  smoothing  and 
surface  reconstruction.  This  process  controls  both  mesh  stretching  and  cell  skewing  while 
maintaining  an  accurate  representation  of  the  surface  geometry. 

3.1  Detection/Division 

The  feature  detection  algorithm  presently  used  is  essentially  a  refinement  of  the  technique  used  in 
the  preliminary  2D  investigations  in  Ref.  [2J.  This  algorithm  was  extended  to  3D  viscous  flow 
in  Ref.  [15]  and  is  fully  documented  in  Ref.  [19],  The  process  consists  of  several  steps  which 
intend  to  isolate  regions  of  locally  high  truncation  error  and  then  reduce  this  error  through  cell 
division.  After  first  examining  the  domain  with  an  undivided  second  difference  of  pressure  to 
find  “shock”  cells,  the  routine  re-scans  only  the  remaining  cells  in  the  computation  to  locate 
“smooth”  features.  This  second  scan  relies  upon  an  undivided  second  difference  of  velocity 
magnitude  to  locate  regions  of  rapidly  varying  shear  stress  and  an  undivided  first  difference  of 
density  to  locate  inviscid  features  in  the  flow.  The  cells  tagged  in  the  second  scan  are  evaluated 
on  a  statistical  basis  as  outlined  in  Ref.  [23].  The  division  routine  operates  recursively  in  a 
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direction-by-direction  manner  and  retains  the  ability  to  divide  the  ceils  directionally  when  features 
are  nearly  mesh  aligned. 

3.2  Mesh  Smoothing 

After  the  division  process  is  complete,  the  entire  mesh  is  smoothed  through  the  sequential 
application  of  Laplacian  and  vector  cross  operators  to  reduce  local  stretching  and  skewing  errors 
within  the  mesh.  This  step  is  important  since  the  adaptation  process  will  reduce  the  truncation 
error  in  the  solution  only  if  the  mesh  remains  smooth  as  it  refines. 

Since  they  involve  only  first  derivative  terms,  discrete  solutions  of  purely  inviscid  flow  problems 
are  more  tolerant  of  poor  meshes.  In  2D,  for  example,  using  trapezoidal  integration  to  evaluate 
the  inviscid  flux  integral  in  Eqs  .{8}  and  {9}  one  may  guarantee  at  least  first-order  accurate 
nodal  updates  (consistency)  regardless  of  mesh  quality^28!.  Viscous  simulations,  however, 
require  computing  the  first  derivative  of  first  derivatives  -  a  process  which  compounds  numeri¬ 
cal  errors  due  to  mesh  stretching  or  skewing.  This  fact  makes  it  impossible  to  guarantee  a 
consistent  representation  of  the  Navier-Stokes  equations  on  an  arbitrary  mesh  (assuming  linear 
basis  functions).  These  observations  lead  to  the  usual  requirements  of  smoothly  varying,  non- 
skewed  computational  cells  for  accurate  representation  of  viscous  flows. 

The  elliptic  smoothing  process  is  separated  into  two  relaxation  operations.  First,  a  Laplacian 
operator  sweeps  through  the  mesh  nodes,  relaxing  a  given  node  toward  the  mid-point  location 
defined  by  the  geometric  center  of  its  neighbors.  This  sweep  proceeds  on  a  direction-by¬ 
direction  basis.  The  operator  degenerates  naturally  on  boundary  faces  or  edges,  so  that  such 
nodes  remain  free  to  slide  along  the  face  or  edge. 

After  each  pass  with  this  smoother,  the  routine  applies  a  vector  cross  operator  to  each  physical 
cell  in  the  mesh.  For  each  mesh  cell,  this  operator  computes  a  right  parallelepiped  with  the  same 
median  dimensions  and  uses  a  relaxation  procedure  to  drive  the  comers  of  the  physical  cell 
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toward  those  of  the  right  parallelepiped.  In  this  manner,  cell  skewing  is  reduced  throughout  the 
mesh. 

Although  nodes  are  free  to  slide  on  and  off  boundary  surfaces,  neither  smoother  alters  the  mesh 
connectivity.  After  each  smoothing  pass,  the  algorithm  checks  to  ensure  that  all  cell  volumes  are 
physically  meaningful.  The  smoother  operates  in  about  10%  of  the  time  required  for  one 
solution  iteration,  and  uses  typical  relaxation  factors  of  approximately  0. 1  with  10  steps. 

3.3  Geometry  Definition 

Since  the  smoother  allows  mesh  nodes  to  move  along  boundaries,  surface  nodes  must  be 
constrained  to  the  body  surface.  This  operation  requires  interrogation  of  a  surface  geometry  data 
base.  Such  a  data  base  also  ensures  that  new  points  introduced  on  the  surface  through  adaptation 
are  placed  directly  upon  that  boundary. 

Figure  5  presents  a  situation  common  to  all  adaptive  schemes  regardless  of  control  volume 
shape.  As  shown  in  the  figure,  when  a  new  point  is  introduced  on  a  mesh  boundary,  it  may  not 
fall  directly  upon  the  body  surface.  The  high  aspect  ratio  cells  typically  found  in  Navier-Stokes 
simulations  exacerbate  this  problem,  and  near  a  highly  convex  surface,  several  new  nodes  may 
be  introduced  which  are  actually  inside  the  surface.  To  rectify  this  situation,  all  surface  nodes 
are  moved  to  the  actual  surface  and  a  search  algorithm  then  follows  each  computational 
coordinate  (for  as  long  as  it  exists)  into  the  mesh  to  adjust  interior  nodes  which  may  have  been 
wrongly  introduced  inside  the  body. 

Several  approaches  exist  to  relocate  the  surface  points^29!.  The  present  algorithm  adopts  an 
approach  which  requires  only  an  unconnected  surface  point  distribution  on  the  body.  This  data 
base  may  be  structured  or  unstructured  and  processing  begins  with  a  triangular  tessellation  of  the 
surface  data  to  form  triangular  patches  describing  the  surface  geometry.  These  patches  are 
currently  planar,  although  no  fundamental  difficulty  exists  in  using  a  higher  order  description. 
Surface  reconstruction  begins  by  projecting  a  ray  from  each  boundary  node  toward  the  surface 
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Before  adaptation  After  adaptation  After  adaptation  and 

surface  reconstruction 


Figure  5.  Reconstruction  of  surface  geometry  from  surface  database. 


21 


mesh.  The  algorithm  then  finds  the  triangle  which  is  pierced  by  the  projection  ray  and  relocates 
the  surface  node  to  the  point  of  intersection  of  the  patch  and  the  ray.  Surface  reconstruction  is 
performed  during  alternating  sweeps  of  the  entire  mesh  smoother.  Figure  6  demonstrates  the 
process  using  a  hemispherical  surface,  and  an  initially  cubic  arrangement  of  physical  cells.  The 
example  shows  results  after  only  one  iteration  of  the  complete  mesh  smoothing  algorithm. 
Successful  experiments  have  mapped  meshes  onto  a  variety  of  convex  and  concave  bodies 
including  the  examples  in  the  next  section  and  the  3D  cropped  delta  wing  in  Refs.  [15]  and  [19]. 
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Figure  6.  Triangulated  surface  geometry  and  resulting  computational  mesh  after  mapping 
physical  cells  onto  surface.  Original  mesh  consisted  of  14x14  cube.  1  smoothing  pass,  and  1 
surface  reconstruction  sweep. 
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4.  Numerical  Experiments  and 

Discussion 


The  numerical  results  contained  within  this  section  intend  to  focus  the  discussion  on  issues 
raised  in  the  introduction  and  to  provide  an  account  of  the  scheme's  behavior.  Herein,  attention 
is  restricted  to  two-dimensional  example  problems.  Refs.  [15]  and  [19]  present  three 
dimensional  results  with  the  current  method  including  comparison  with  experiment 

4.1  AGARD  03  Test  Case 

Before  examining  the  behavior  of  the  full  viscous  scheme,  it  is  appropriate  to  establish  the 
accuracy  of  the  convective  discretization  and  adaptive  methodology.  The  AGARD  03  test  case 
has  been  explored  in  the  literature  with  a  variety  of  methodsn2],[30],[31]  it  examines  symmetric, 
inviscid,  flow  at  Mach  0.95  over  a  NACA  0012  profile  and  the  case  presents  a  particular 
challenge  to  numerical  schemes!12!  since  it  contains  strong  coupling  between  smooth  and 
discontinuous  features  which  span  many  length  scales.  The  flow  becomes  supersonic  over  the 
airfoil  and  forms  an  oblique  ‘’fish-tail”  shock  system  attached  to  the  trailing  edge.  A  weak 
normal  shock  forms  downstream  of  these  shocks,  in  the  wake  of  the  airfoil.  The  smooth  ex¬ 
pansion  over  the  airfoil  surface  weakens  and  bends  the  oblique  shocks  away  from  the  body. 
Finally  a  shock  triple  point  forms  which  locates  the  normal  shock  in  the  wake. 

The  shock  polar  at  the  7.99°  trailing  edge  of  the  airfoil  is  extremely  steep  at  these  Mach  numbers 
and  small  changes  in  the  trailing  edge  Mach  number  may  radically  alter  the  fish-tail  shock 
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structure!  121.  This  sensitivity  is  further  compounded  by  the  fact  that  the  shock  triple  point  lies 
approximately  5  chords  away  from  the  airfoil.  Stiffness  in  the  crossflow  direction  makes  this 
case  quite  sensitive  to  far  field  boundary  placement.  Successful  prior  calculations  located  the 
boundary  about  100  chords  away!12!*!31!.  Furthermore,  fine  resolution  of  the  shocks  in  the 
wake  does  not  ensure  low  overall  truncation  error  since  small  errors  in  smooth  regions  of  the 
flow  may  cause  gross  errors  in  shock  location,  regardless  of  resolution!  121.  For  these  reasons, 
the  problem  directly  addresses  the  issues  of  mesh  convergence  and  consistency,  raised  earlier, 
while  simultaneously  evaluating  the  inviscid  discretization. 

The  Mach  contours  in  Figure  7  provide  an  overview  of  the  discrete  solution.  In  this  simulation 
the  entropy  cutoff  in  the  TVD  scheme  was  applied  to  all  eigenvalues.  The  figure  shows  both 
Mach  contours  and  the  final  adapted  half-mesh  with  5  levels  of  mesh  cells  and  approximately 
18,000  nodes.  The  far  field  boundary  in  the  calculation  was  located  75  chords  away,  and  used 
the  characteristic  based  treatment  of  Ref.  [32].  Figure  8  shows  an  enlargement  of  the  Mach 
contours  in  the  region  near  the  airfoil. 

Mesh  refinement  studies  in  Ref.  [12]  using  CFL2D  and  the  floating  shock  Fitting  procedure  of 
Ref.  [31]  provide  a  basis  for  evaluating  the  discrete  solution.  The  finest  mesh  in  these  studies 
consisted  of  1,567,485  nodes  (2049  x  765)  and  the  results  from  these  studies  located  the  normal 
shock  (measured  to  the  sonic  line)  between  3.32c  and  3.35c  downstream  of  the  trailing  edge. 
The  current  results  place  this  structure  at  3.34±0.02c.  The  Mach  number  of  the  flow  through  the 
oblique  shocks  varied  from  1.45  near  the  surface  to  1.25  near  the  shock  triple  point.  These 
values  also  match  those  from  CFL2D. 

4.2  Flat  Plate  Boundary  Layer  Flow 

The  viscous  test  cases  begin  with  a  simulation  of  flow  over  a  flat  plate  at  Mach  0.5  and  Reynolds 
number  of  5,000  per  unit  length.  The  case  was  run  on  unadapted  meshes  extending  2.5 L  above 
and  in  front  of  the  plate.  Ref.  [19]  considers  a  similar  example  using  the  current  method. 
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Figure  7.  Mach  contours  (inc.=0.05)  and  adapted  upper  half  mesh  for  AGARD  03  test  case. 
Mach  0.95,  a  =  0°,  inviscid  flow.  18,000  nodes. 


Figure  8.  Near  field  Mach  contours  (inc.=0.05)  for  AGARD  03  test  case. 
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However,  the  present  discussions  warrant  reexamination  of  this  case  to  provide  insight  into  the 
scheme’s  behavior  when  the  entropy  cutoff  is  applied  to  all  the  eigenvalues. 

The  Upwind  TVD  convective  formulation  makes  use  of  an  entropy  cutoff  which  precludes 
violation  of  the  entropy  condition  at  points  where  the  eigenvalues  vanish  in  the  flown7],[2].[l9] 
Since  all  linear  eigenvalues  become  zero  at  solid  walls,  this  cutoff  may  adversely  affect  the 
scheme's  accuracy  in  such  regions.  Perhaps  the  most  common  approach  to  avoiding  this 
situation  is  to  apply  the  cutoff  only  to  the  nonlinear  eigenvalues.  Figure  9  presents  u-velocity 
profiles  plotted  versus  similarity  parameter  and  compared  directly  with  a  Blasius  profile.  Two 
cases  are  examined,  one  in  which  all  the  eigenvalues  were  thresholded,  and  a  second  one  which 
applied  the  cutoff  only  to  the  nonlinear  waves.  Since  these  effects  are  most  obvious  with 
relatively  sparse  boundary  layer  resolution,  this  example  uses  only  about  six  points  in  the  layer. 
Despite  this  apparent  lack  of  resolution,  when  the  cutoff,  5,  applied  only  to  the  nonlinear  waves, 
the  solution  already  provides  a  reasonable  approximation  of  the  Blasius  profile.  These  results 
are  typical  of  those  obtained  with  high-resolution  upwind  methods  and  often  lead  to  the  practice 
of  cutting  off  only  the  nonlinear  waves  in  viscous  flow  simulations^33!-!34!. 

Previous  results  with  the  present  scheme!1 51*U9]  follow  this  approach  and  apply  S  only  to  the 
nonlinear  waves.  Figure  10,  however,  presents  u  and  v-profiles  obtained  cutting  all  the  waves 
in  the  solution,  and  using  cutoff  values  of  8=  0.01,  0.1,  and  0.3.  By  increasing  the  boundary 
layer  resolution  to  about  13  points,  the  detrimental  affects  of  Jon  the  strcauiwise  velocity  nearly 
vanishes.  For  comparison.  Figure  1 1  contains  the  same  u  and  v-velocity  profiles  computed  with 
the  central  difference  option  and  4th  difference  artificial  dissipation  (V2  =  0.0,  v*=  1/128  - 
1/32). 

With  13  points  in  the  boundary  layer,  agreement  with  the  Blasius  profile  is  excellent  with  both 
the  central  and  TVD  discretization  of  the  inviscid  terms.  However,  close  examination  reveals 
some  slight  discrepancies.  These  appear  predominantly  in  the  central  difference  results  at  V4  = 
I/32  and  */64-  °f  these  curves  display  some  evidence  of  a  "viscous  overshoot”  just  outside 
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Figure  9.  Effects  of  cutting  off  all  or  only  the  nonlinear  eigenvalues  on  flat  plate  boundary  layer 
simulation  using  the  TVD  discretization.  Comparison  with  Blasius  similarity  solution,  = 
0.5,  ReL-  5,000. 
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Figure  10.  Velocity  profiles  for  flat  plate  boundary  layer  solutions  using  TVD  convective 
discretization  and  applying  the  entropy  cutoff  to  all  eigenvalues.  Above:  u-velocity  profiles. 
Below:  v- velocity  profiles,  using  3  different  values  for  the  cutoff  8. 
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Figure  11.  Velocity  profiles  for  flat  plate  boundary  layer  solutions  using  central  difference 
convective  discretization  with  4th  difference  artificial  dissipation  {v4).  Above:  u-velocity 
profiles.  Below:  v- velocity  profiles,  using  v4  =  l/32,  v4  =  l/64,  and  v4  = 
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the  knee  of  the  profile.  The  excessive  dissipation  in  these  simulations  does  not  permit  the  profile 
to  bend  rapidly  enough  to  follow  the  Blasius  profile  around  this  knee.  This  error  decreases 
rapidly  with  decreasing  V4,  and  almost  no  evidence  of  it  is  present  in  the  simulation  at  V4  = 
Vl28-  The  u-velocity  profiles  of  the  TVD  discretization  show  almost  no  variation  with  the  value 
of  the  cutoff,  S.  Additionally  there  is  no  evidence  of  the  viscous  overshoot,  simply  because  the 
flux  limiter  would  see  such  an  overshoot  as  a  “local  maximum”  and  clip  it  off. 

Slightly  more  variation  is  seen  between  the  solutions  in  v-  velocity  profiles  in  Figs.  10  and  1 1. 
Here  the  viscous  overshoot  in  the  central  schemes  is  much  more  readily  apparent.  Although  the 
effects  of  the  entropy  cutoff  are  somewhat  more  apparent  in  these  curves,  the  profiles  for  the 
TVD  still  remain  quite  similar  to  one  another . 

Figure  12  shows  skin  friction  development  along  the  flat  plate.  This  plot  reports  results  for  both 
the  central  and  TVD  discretizations  at  two  levels  of  resolution  and  compares  these  with  the 
Blasius  relation  C/=  0.664(/?ex)1/2.  With  13  points  in  the  boundary  layer  (at  a  local  Rex  of 
5000),  the  theoretical  skin  friction  is  predicted  well  by  about  Rex  =  100,  or  2%  of  the  plate 
length  from  the  leading  edge.  With  26  points  in  the  layer,  the  skin  friction  matches  almost 
immediately.  With  all  the  eigenvalues  cutoff,  lower  resolution  yields  generally  poor  agreement 
between  the  TVD  and  Blasius  result. 

4.3  Symmetric  Laminar  Airfoil  Test  Case 

The  final  example  considers  symmetric  viscous  flow  over  a  NACA  0012  at  a  Mach  number  of 
0.5  and  Reynolds  number  of  5,000  with  adiabatic  wall  boundary  conditions.  The  recent 
literature  includes  investigations  of  this  flow  by  a  host  of  structured  and  unstructured 
methodsl35!’!36!  and  the  case  is  well  understood.  In  addition  to  these  computational  results,  the 
relative  simplicity  of  the  flow  permits  an  analytic  analysis  which  provides  an  independent  check 
of  self-consistency  in  the  discrete  solutions. 
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As  it  is  given  in  Ref.  [37],  the  airfoil  profile  does  not  quite  close  at  the  trailing  edge.  This 
shortcoming  was  circumvented  by  continuing  the  profile  until  reaching  the  centerline  ard  then 
rescaling  both  axes  to  place  the  trailing  edge  again  at  1.0.  Such  detail  is  important  in  this  flow 
since  the  case  is  near  the  upper  boundary  for  stable  laminar  flow.  Additionally,  at  the  conditions 
of  this  simulation,  a  separation  bubble  appears  at  the  trailing  edge,  and  the  separation  location 
depends  directly  on  the  geometry  of  the  trailing  edge. 

Figure  13  displays  a  view  of  the  final  adapted  mesh  in  the  vicinity  of  the  airfoil.  The  starting 
mesh  for  this  calculation  was  a  25x15  half  C-mesh,  and  the  final  adapted  mesh  contained  5  levels 
of  cells  and  8,100  nodes.  Two  inset  regions  on  this  mesh  show  the  finest  two  levels  of 
adaptation  near  the  leading  and  trailing  edges. 

One  particular  challenge  presented  by  this  test  case  lies  in  properly  resolving  the  extent  of  the  thin 
separation  bubble  at  the  trailing  edge.  In  resolving  this  structure,  the  adaptation  routine  has 
changed  the  cell  aspect  ratio  from  approximately  20: 1  to  only  5: 1  near  the  separation  point.  The 
final  two  adaptations  divided  cells  primarily  in  the  streamwise  directions.  As  a  result  the  cell 
aspect  ratio  on  the  wall  varies  from  about  50  at  the  midchord  to  only  5  near  the  separation 
location.  It  is  interesting  to  note  that  in  the  nonadaptive  efforts  in  Refs.  [35]  and  [36]  the  cells 
had  approximately  the  same  streamwise  spacing,  but  used  between  1/2  and  1/5  the  wall  normal 
mesh  spacing.  The  current  procedure  arrived  at  this  streamwise  spacing  by  adaptation  without 
user  intervention. 

Figure  14  presents  contours  of  constant  Mach  number  near  the  airfoil  and  an  enlargement 
showing  computed  streamlines  near  the  trailing  edge.  Figure  15  displays  the  variation  in  surface 
pressure  and  skin  friction  coefficients  over  the  airfoil.  These  profiles  agree  well  with  established 
previously  published  resultsf35TC3<5J_ 
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Figure  12.  Sirin  friction  evolution  along  flat  plate  using  both  Upwind  TVD  and  central 
differencing.  Entropy  cutoff  applied  to  all  eigenvalues  in  TVD  solutions.  Blasius  solution 
provided  for  comparison.  Af„  =  0.5,  Rei  =  5,000. 
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Figure  13.  Adapted  half-C  mesh  near  airfoil  for  viscous  NACA  0012  test  case.  5  levels  of 
cells,  8,100  nodes,  =  0.5,  Rec  =  5,000. 
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Figure  14.  Mach  contours  (inc.  =  0.01)  in  discrete  solution  using  TVD  discretization  and 
applying  entropy  cutoff  to  all  eigenvalues.  Separation  point  at  8 1 .5%  chord.  =  0._  ,  Rec 
=  5,000. 
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4.4  Separation  Location 

Table  1  compares  the  location  of  the  separation  point  as  computed  with  several  methods.  Figure 
16  shows  streamlines  and  velocity  vectors  near  the  separation  region,  and  has  the  half  saddle  of 
separation  labeled  on  the  surface.  As  indicated  in  the  table.  Ref.  [35}  showed  that  this  location 
was  sensitive  to  the  level  of  numerical  dissipation  in  the  solution. 

With  the  current  method,  dissipation  enters  from  the  entropy  cutoff  and  from  the  discretization 
error.  The  adaptation  present  in  the  final  solution  aims  to  reduce  this  error  by  refining  the  mesh. 
Since  the  separation  location  is  sensitive  to  the  level  of  numerical  dissipation,  it  is  useful  to 
examine  its  behavior  as  the  mesh  refines.  Figure  17  traces  the  separation  location  on  the  last  four 
meshes  used  in  the  solution  process.  As  the  background  dissipation  is  reduced  by  mesh  re¬ 
finement,  the  separation  point  moves  to  a  location  consistent  with  a  solution  of  the  governing 
equations  (since  the  modified  PDE's  approach  the  true  governing  equations).  It  is  interesting  to 
see  the  extent  to  which  the  last  two  adaptations  affected  the  location  of  the  separation  point,  since 
the  cells  were  divided  only  in  the  streamwise  directions  and  the  wall  normal  mesh  spacing 
remained  constant. 


Although  the  curve  in  Fig.  17  appears  to  flatten,  it  does  not  actually  asymptote  convincingly  to 
the  value  of  81.5%  reported  earlier.  Equation  {18}  provides  a  means  to  further  investigate  the 
magnitude  of  the  discretization  error  in  the  immediate  vicinity  of  the  separation  point. 


,,  deo  / 

tan(0)=  ^  dx  dp 

/  dx 


{18} 


Table  1.  Separation  locations  predicted  by  several  methods. 


Method 

Dissipation 

Sep.  Loc. 
<%d 

Wall  spacing 

Present  Method 

(5  =  0.1 

81.5 

0.001 

Triangular  scheme  of  [35] 

v<=l/256 

81.4 

0.0002 

Triangular  scheme  of  [35] 

v*= 1/64 

83.4 

0.0002 

Central  Diff  Method  of  [26] 

N/A 

81.9 

N/A 

Central  Orff  Method  of  [361 

N/A 

81.1 

0.0006 
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This  expression  for  the  classical  two-dimensional,  laminar  separation  angle,  6,  was  derived  in 
Ref.  [38],  In  Ref.  [39]  the  authors  arrive  at  the  same  expression  from  critical  point  theory  where 
they  show  that  the  slope  of  the  isocline  emanating  from  a  half-saddle  of  separation  may  be 
derived  from  a  direct  balance  between  the  normal  and  tangential  momentum  equations.  This 
facilitates  a  direct  comparison  between  the  measured  angle  in  Fig.  16  and  the  angle  predicted  by  a 
local  solution  of  the  (incompressible)  Navier-Stokes  equations.  The  measured  separation  angle 
in  Figure  16  is  6.07±0.2°,  while  evaluating  Hq.  {18}  numerically  results  in  0  =  6.09±0.1 4°. 

This  analysis  does  not  verify  the  location  of  the  separation  point,  but  it  does  indicate  that  the 
discretization  error  in  the  region  of  the  separation  is  very  small.  Thus,  the  extent  of  the 
separation  bubble  is  not  “held  back”  by  locally  high,  nonphysical,  dissipation  or  other  dis¬ 
cretization  error.  For  comparison,  performing  the  same  analysis  on  the  coarsest  mesh  in  Fig.  17 
(sep.  location,  at  96%)  yields  a  9.5°  discrepancy  between  the  measured  and  computed  isocline 
angles,  indicating  large  discretization  errors. 

4.5  Entropy  Cutoff  and  Non-Physical  Solutions 

The  preliminary  computation  of  this  viscous  NACA  0012  test  case  adopted  the  usual  approach  of 
applying  the  entropy  cutoff  to  only  the  nonlinear  waves  in  viscous  simulate.  s.  While 
converging  after  the  second  adaptation  (3rd  mesh),  the  method  resulted  in  solutions  similar  to 
that  shown  at  the  right  of  Figure  18.  Such  nonphysical  solutions  persisted  for  CFL  numbers 
from  0.1  to  4.3  and  for  3  and  5  stage  Runge-Kutta  schemes,  using  local  or  global  time  steps. 
The  problems  were  considered  anomalous  since  the  same  code  has  performed  well  in  a  variety  of 
two-  and  three-dimensional  test  cases[l5hl19].  The  problem  was  solved  by  applying  the  entropy 
cutoff  to  all  wave  speeds,  and  it  is  this  experience  which  motivated  the  inclusion  of  the  flat  plate 
investigation  with  all  the  eigenvalues  cut 

The  appearance  of  such  behavior  in  an  otherwise  well  behaved  and  well-documented  solver 
warrants  some  discussion.  Similar  anomalous  results  have  been  documented  by  several  authors 
133], [34], (40J, [41]  Ruling  out  bugs  in  the  code  (based  upon  previous  success  with  a  wide  variety 
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Figure  18.  Stagnation  region  at  leading  edge  of  NACA  0012  after  second  adaptation.  Frame  a 
entropy  cutoff  of  8=  0.1  applied  to  all  eigenvalues.  Frame  b:  entropy  cutoff  of  8-  0. 1 
applied  to  nonlinear  eigenvalues  only. 
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of  test  cases),  the  two  most  likely  explanations  are:  1)  problems  in  the  TVD  scheme  itself,  2) 
problems  resulting  from  the  interaction  of  the  TVD-inviscid  and  central-viscous  discretizations. 

Close  examination  of  Fig.  18b  shows  that  the  boundary  layer  near  the  stagnation  point  remains 
well  behaved.  It  is  only  in  the  inviscid  stagnation  region  outside  this  viscous  layer  that  problems 
arise.  In  that  region,  the  viscous  terms  are  smaller  by  2-3  orders  of  magnitude  than  near  the 
“knee”  of  the  boundary  layer  profile.  This  statement  favors  explanation  #1,  and  shifts  focus  to 
the  TVD  scheme  itself.  The  only  difference  in  the  inviscid  discretization  used  in  Fig.  18b  and  a 
standard  inviscid  flow  (like  the  AGARD  03  test  case)  is  the  fact  that,  in  18b,  the  linear 
eigenvalues  in  the  viscous  flow  were  not  thresholded  near  zero. 

The  cutoff,  S,  is  introduced  into  the  method  to  prevent  violation  of  the  entropy  condition!16!. 
When  the  cutoff  is  invoked,  it  introduces  (2nd  difference  like)  explicitly  added  dissipation.  The 
added  dissipation  term  increases  the  entropy,  and  only  physically  consistent  results  should 
remain  in  the  space  of  possible  discrete  solutions.  In  the  stagnation  region,  all  the  linear  eigen¬ 
values  go  to  zero.  In  the  inviscid  part  of  this  stagnation  region,  if  no  dissipation  is  introduced, 
the  preference  for  selecting  a  solution  consistent  with  the  entropy  condition  is  lost.  On  a  coarse 
enough  mesh,  the  dissipation  introduced  through  the  discretization  error  may  be  sufficient  to  rule 
out  such  inconsistent  discrete  solutions.  However,  when  the  mesh  is  refined  (as  by  adaptation), 
such  dissipation  will  diminish  at  a  rate  determined  by  the  order  of  accuracy  of  the  scheme. 
Therefore,  on  a  fine  enough  mesh,  nonphysical  solutions  may  appear.  Such  a  hypothesis 
explains  the  converged  and  physically  consistent  results  obtained  on  the  first  two  meshes  in  the 
adaptation  sequence.  Notice  that  inside  the  boundary  layer,  u  also  tends  to  zero.  Here, 
however,  the  viscous  terms  in  the  Navier-Stokes  equations  generate  entropy  with  the  apparent 
result  of  “selecting”  physical  solutions.  This  statement  is  supported  by  the  seemingly  well 
behaved  velocity  vectors  seen  just  inside  the  boundary  layer  in  Fig.  18b. 

Further  supporting  evidence  for  this  hypothesis  came  from  examining  the  region  near  the  saddle 
point  of  separation.  Near  the  point  of  inflection  in  the  velocity  profile,  the  virtual  lack  of 
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curvature  gives  rise  to  generally  small  viscous  stresses.  When  this  point  is  near  the  wall,  the 
velocities  are  also  small  and  so  the  eigenvalues  will  tend  to  vanish  and  discretization  may  violate 
the  entropy  condition.  Investigation  of  the  flow  topology  in  this  area  revealed  physically 
inconsistent  streamline  patterns  in  the  vicinity  of  the  separation  point  just  away  from  the  wall. 
These  arguments  are  not  intended  as  proof,  but  do  offer  avenues  for  further  investigation. 

This  problem  in  inviscid  stagnation  regions  reemphasizes  the  need  for  thresholding  of  the  wave 
speeds  in  schemes  which  may  violate  the  entropy  condition.  In  complex  flows,  errors  such  as 
the  inconsistent  flow  topology  cited  here  could  easily  go  unnoticed  and  give  rise  to  unpredictable 
flow  physics.  This  statement  goes  against  the  common  practice  of  thresholding  only  some  of  the 
waves,  or  different  waves  in  different  directions  in  the  domain^], [40], [19]  While  such 
techniques  may  produce  excellent  results  for  skin  friction  and  heat  transfer  they  still  do  not,  in 
general,  ensure  convergence  to  only  physically  consistent  solutions. 
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5.  Summary  and  Conclusions 

This  work  presented  further  details  of  a  newly  developed  unstructured  discretization  for 
modeling  the  viscous  terms  in  the  Navier- Stokes  equations.  A  detailed  description  of  this 
discretization  and  its  implementation  was  presented  *n  a  form  requiring  only  cell-to-node  and 
node-to-cell  communication.  The  method  yields  second-order  accuracy  on  smooth  meshes,  and 
is  free  stream  preserving  on  any  arbitrary  mesh.  Various  aspects  of  the  technique  were  evaluated 
through  test  problems,  analysis  and  mathematical  proof.  The  method  was  incorporated  into  a 
previously  documented  adaptive  TVD  scheme  using  hexahedral  based  unstructured  meshes 
which  adapt  to  the  evolving  solution  through  directional  cell  division.  The  discussion  also 
outlined  the  techniques  developed  for  controlling  mesh  quality  and  reconstructing  the  correct 
surface  geometry  after  adaptation  and  mesh  smoothing. 

The  issue  of  mesh  convergence  was  investigated  through  examination  of  the  AGARD  03  test 
case  which  is  known  to  be  particularly  challenging  to  adaptive  schemes.  The  new  scheme 
matched  results  from  previously  reported  mesh  convergence  studies  conducted  with  two 
structured  schemes.  These  discrete  solutions  also  agreed  with  the  results  of  two  adaptive 
upwind  schemes  which  used  a  feature  detection  algorithm  specifically  designed  in  response  to 
the  subtlety  of  the  AGARD  03  test  case. 

A  detailed  examination  of  viscous  flow  over  a  NACA  0012  at  and  Rec  =  5,000  resulted 

in  Cp  and  C/distributions  which  agreed  well  with  previously  published  solutions  on  unadapted 
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meshes.  In  addition,  the  separation  location  was  predicted  in  agreement  with  previous  findings 
and  further  supported  by  an  analysis  of  the  theoretical  laminar  separation  angle. 

During  the  course  of  the  investigations,  it  became  evident  that  the  common  practice  of 
thresholding  only  the  nonlinear  eigenvalues  when  using  upwind  methods  whose  formulation 
permits  violations  of  the  entropy  condition  can  admit  nonphysical  discrete  solutions.  Such 
violations  are  likely  in  stagnation  regions  and  near  points  of  inflection  in  the  boundary  layer 
velocity  profile.  Such  solutions  were  examined  and  related  to  similar  anomalous  solutions 
reported  by  a  variety  of  other  researchers  also  using  Roe  type  schemes.  Since  many  of  these 
nonphysical  solutions  may  go  undetected  in  complex  flow,  this  research  advocates  the 
(admittedly  conservative)  position  of  always  thresholding  both  the  linear  and  nonlinear 
eigenvalues  in  all  directions.  While  this  approach  does  slightly  increase  the  resolution 
requirements  for  accurate  representation  of  the  boundary  layer,  it  appears  to  be  the  only  sure  way 
to  avoid  admitting  nonphysical  solutions  into  the  solution  space.  A  detailed  evaluation  of 
discrete  solutions  to  a  flat  plate  boundary  layer  was  included  showing  good  skin  friction 
prediction  with  approximately  13  points  despite  thresholding  all  of  the  eigenvalues. 
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Appendix: 

Preservation  of  Uniform  Flow  on  an 
Arbitrary  Mesh 

Claim:  AUvisc  remains  identically  zero  on  any  arbitrary,  nonoverlapping  hexahedral 

mesh  when  no  spatial  gradients  arc  present  in  the  state  vector. 

Proof:  Let  G  be  a  3D  spatial  graph  with  a  hexahedral  based  tessellation  containing  node 

i.  Taking  the  surface  vectors  of  the  polyhedra  as  described  in  the  text,  the  viscous  update  to 
node  i  becomes  (Eq.  {11}): 

{Uvuc)i  =  *SN-  Fv'  Ss+  Fv‘Se-  Fv  *  Sw 

V  i 

+  Fv*Sf-Fv-Sb]  (a  .i) 

For  the  Navier-Stokes  equations  under  the  assumptions  stated  in  the  text,  (Fv)M.s.e.w.f.b  are 
linear  combinations  of  the  First  derivative  quantities  dfa  (Eq.{  12}).  Thus,  when  all  =  0 

then  the  viscous  fluxes  (Fv)nse  wfb  =  0  (see  Eqs.{3}-{7})  and  A Uvisc  will  vanish  by 
Eq.{A.  1 }. 
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Referring  to  Figs.  2,  3,  and  4,  we  begin  by  examining  the  members  of  the  tensor  along 
edge  Tj  which  passes  through  the  E  face  of  the  auxiliary  ceil  surrounding  node  i. 

( ds<bi)i  =  2V^+  ~  +  ss?) 

+  ifistf+Sfl-iffaSi  +  Stf)  (A. 3} 

+  0q(Szf  +  s£j)~  <$e(s^  +  s&  )] 

In  a  zero  gradient  field,  however,  are  constant  at  their  reference  (°°)  values. 

(  ~  2^[+  (5?/+  S«D“ 

'  +(^+^)-(s4f+Sj)  {A.4} 

+  (^f +  +  ^f)] 

Thus,  for  the  elements  of  d^v  to  be  zero,  the  bracketed  term  in  {A.4}  must  sum  to  zero.  This 
expression  is  exactly  the  sum  of  the  surface  vectors  of  the  secondary  cells.  Rearranging  the 
bracket  in  this  expression,  we  seek  to  evaluate: 

l+Stf+Sj-Sf-S'r+Sfi+Sjl 
-  sg-s£  +  S?f  +  s£  -  sg  -  Sjj  ] 

Collecting  terms  yields: 

[+sjstr+sfi-st  +  s£-sj 

+SSJ  ~S^  *  StJ~  S(J  +  S(J  -  Stl  i  { A.5 ) 

but  Sq  =  Sfy  because  it  is  the  common  face  between  auxiliary  cells  /  and  j  and  the  expression 
in  {A.5}  may  be  re-written  as: 

[+  S^i  -S^i  +  Sy  -  Sfr  +  S^i  -  %  { A .6} 

+S$j~SsJ  +  S£~  S£  +  Sfj  ~  S£  ] 

We  note  that  when  the  auxiliary  cells  are  constructed  as  described  in  the  text,  they  form  a  dual 
mesh  consisting  of  closed  volumes  which  fill  the  interior  space  of  the  graph  without  overlapping 
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or  interior  voids.  Since  the  auxiliary  cells  are  all  closed  volumes,  their  surface  vectors  sum  to 
zero,  and  we  have: 


S£  +  s£  ~  s£  +  S£  ~s£  =°  { A. 7 } 

sf-sf  +  sj-sj  +  sj-sj.  o 

Expression  { A.6}  then  becomes  identically  zero,  which  ensures  that  the  sum  of  surface  vectors 
in  Eq.{A.4}  will  always  be  zero.  Thus,  in  accordance  with  the  statements  following  Eq.{A.l }. 
AUVisc  will  always  be  identically  zero  in  the  presence  of  an  unperturbed  free  stream  flow. 
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